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Abstract 

We study exact results concerning the non-affine displacement fields observed by Tanguy et al 
[Europhys. Lett. 57, 423 (2002), Phys. Rev. B 66, 174205 (2002)] and their contributions to 
elasticity. A normal mode analysis permits us to estimate the dominant contributions to the non- 
affine corrections to elasticity and relate these corrections to the correlator of a fluctuating force 
field. We extend this analysis to the visco-elastic dynamical response of the system. 
Keywords: amorphous solids, Born-Huang approximation, visco-elasticity, non-affine 
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A straightforward estimate of the elastic constants of simple crystals can be performed 
in the (classical) zero temperature limit: the relative initial positions of atoms are known; 
elementary deformations are homogeneous even at the microscopic level. It is thus a simple 



task to add up all contributing interactions. These assumptions-zero temperature anc i 
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mogeneous displacement of the particles-constitute the basis of the Born-Huang theory. 
These assumptions can also be used to estimate the elastic constants of a disordered struc- 
ture: they provide approximate expressions involving integrals over the pair correlation; in 
liquid theory, these expressions correspond to the infinite frequency moduli. Of course, 
the two assumptions of zero temperature and homogeneous displacement are not valid in 
general and corrections to the Born-Huang approximation are expected to arise from the 
failure of either. 

Early studies by Squire, Holt and Hoover, P, |(| focused on thermal contributions to 
elasticity in crystals. More recently, a surge of interest for athermal materials, like gran- 
ular materials or foams, attracted some attention to corrections to the Born-Huang ap- 
proximation which arise solely from the non-trivial structure of the potential energy land- 
scape. Namely, in disordered solids at zero temperature, the assumption 
that particles follow homogeneous (affine) displacement fields is incorrect: when a mate- 
rial is strained-even by vanishingly small amounts of deformation-particles minimize the 
potential energy of the system by following non-affine displacements. This idea was re- 
cently r ecog nized in numerical simulations of compressed emulsions j3] and Lennard- Jones 
glasses. |lfll. 11 ) In particular, Tanguy et al jlU [l]| have clearly shown that non-affine cor- 



rections to the Born-Huang approximation hold in the continuum limit and amount to an 
important fraction of the Born-Huang term itself. It is thus essential to understand these 
corrections well if we ever want to be able to construct approximations to the elastic con- 
stants of amorphous materials. 

Formal expressions for the non-affine (or "inhomogeneous" , in the language of Wal- 
lace Jl2|) corrections to the Born-Huang approximation were written early on 
[3 llfij- These formal expressions have been used almost exclusively as a tool to calculate 
elastic constants in computer simulations, but were given little attention in light of their 
basic importance. We believe that this arose from two limitations: (i) prior works remained 
at an essentially technical level, aiming merely to provide tools for numerical simulations (ii) 
the derivations of formal analytical expressions for elastic constants have always relied on 
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simplifying assumptions-either dealing with an infinite system or valid at zero stress. These 
simplifications make it difficult to ascertain the domain of validity of various formulae or 
symmetries. In response to these issues, we wish to attack this problem from two opposite 
angles: provide an even more systematic and general treatment than before, yet relate this 
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formalism to numerical observations similar to Tanguy et al. 
treatment could serve, at least, as an introduction to the subject. 

From the formalism, we wish to extract information about the contribution of various 
scales to the non-affine corrections to elasticity: the question we have at heart is whether 
these corrections originate from small or large scales, or involve a broad distribution of scales 
as suggested by the observation of vortex-like patterns. [^] This question is directly relatec. 
to the existence of a continuous limit for elastic properties of amorphous structures, jiol 
To address these questions, we perform a normal mode decomposition of the non-affine 
displacement field: it permits quantifying the contribution of every frequency shell to the 
non-affine corrections to elastic constants. In the large size limit, these contributions seem 
to be self- averaging quantities (in the sense that an ensemble average over subsystems will 
produce results which are equivalent to a single large system). The existence of this self- 
averaging property leads to an expression for the elastic constants in the continuum limit 
which resembles the sum rules of liquid theory. 

Finally, our attention was attracted by several related issues in the recent literature. 
Studies of sound potion aud attenuation in granul a r mater ia ls QQ, HQ 0, of the 
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22j,|23j emphasize the need for 



visco-elastic response of dense emulsions and foams, 
a deep theoretical understanding of the visco-elastic response of disordered solids. Related 
experiments by McKenna and coworkers indicated that features of the visco-elastic response 



of amorphous materials are related to measurable changes in their elastic constants. |24( We 
thus complement the study of static response by a study of dynamical response, and derive 
a formal expression for the visco-elastic moduli. We shall establish that the relaxation 
spectrum is directly and simply related to a correlator emerging from the normal mode 
analysis. 

The present work is meant to present the general framework of our analysis, which will 
be the basis of future numerical studies. Although, here, we will use numerical simulations 
to illustrate analytical developments, the main core of our numerical study will be presented 
in a dedicated article. j3] 
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1. NON-AFFINE DISPLACEMENT FIELDS 



In this work we consider the mechanical response of an amorphous solid quenched at zero 
temperature. Our formalism permits dealing explicitly with finite size systems: it rests on 
the idea that, during a quench at zero temperature, any finite size system relaxes toward one 
of many local minima in the potential energy landscape. j3| Being at zero temperature, the 
system is then prescribed to lie at this minimum at all times. Small external perturbations 
are then expected to induce continuous changes in the local minimum. Large external 
perturbations may induce the vanishing of the local minimum occupied by the system: this 
vanishing occurs when the basin of attraction of this minimum reduces to a single point, 
that is when the minimum collides with at least one saddle point, j^i 28, 2gL Isoj] 



The difference between small, continuous changes of the local configuration and 
catastrophic events is exemplified figure ^ where the response of an amorphous system is 
shown as a function of shear. The parameter 7 measures the total shear deformation from 
a quenched state. In this picture, continuous segments are associated with shear-induced 
changes of local minima and discontinuities to their vanishing. After a minimum has disap- 
peared, the system, coupled to a zero temperature thermal bath, relaxes in search of a new 
minimum in the potential energy landscape. 3(1 [nj 

The separation of continuous changes of the local minimum and catastrophic events is, 
of course, not limited to shear deformation. In a more general context, suppose that we 
denote by 7 the amplitude of any external drive applied to the system. Suppose that the 
system lies at a minimum which, for 7 = 70, is far from any catastrophe. Let us denote 
^7 = 7 — 7o the amplitude of a perturbation of the external drive around 70. If we were to 
slowly increase the amplitude of the perturbation from 0, the system would smoothly follow 
a trajectory 74(7) m configuration space for all amplitudes £7 < e such that the basin of 
the local minimum remains non-vanishing. If the Hessian is non-degenerate, this trajectory 
7^(57) is unique. If we now stop the perturbation at some 57 max < e then revert the drive 
down to £7 = 0, the system would simply follow the trajectory 7^(7) from 7 = 70 + 57max 
to 70, backward. In this sense we will say that the continuous segments are microscopically 
reversible. 

When the system is driven quasi-statically along these continuous segments, no energy is 
dissipated. The reason is that at each point along these segments, the system is at mechanical 
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FIG. 1: A typical energy vs. strain curve obtained in an athermal quasistatic simulation as 
described in the text. Note the smooth elastic segment which terminates at a plastic discontinuity. 

equilibrium: the force applied to any particle is exactly zero and zero forces do no work. 
The quasi-static response corresponds to the thermodynamically reversible, elastic, part of 
the mechanical response. We will see, however, that energy is dissipated when the system is 
driven at finite deformation rates along these continuous segments. This dissipation results 
from the fact that finite deformation rates induce non-zero forces which dissipate energy via 
the coupling with the zero temperature thermal bath. 

1.1 Notation 

We shall let underline and double-underlines respectively indicate vectors and tensors 
referred to a fixed Cartesian system (xi,X2,xs). We shall also use the convention that 



Greek indices refer to Cartesian components of vectors or tensors, while Roman indices refer 
to the particle numbers. Bold type denotes fields which are defined on every particle in the 
material: 

f = {fi}ie{l,...,n} , A = {(A,a^)}ie{l,...,n} • 

Dots and double dots indicate matrix products and summation convention is always applied 
on repeated (Greek and Latin) indices. By convention, we also write: 

fdA\ = (dA a \ aM = ( d 2 A a \ 

\drj af3 \drpj &R \drds) aM \drpdsj 

The superscript T indicates the transpose of a matrix, _1 its inverse, and ~ T the inverse of 

its transpose. The identity matrix is denoted 1, and its components 5 a p. 



1.2 Displacement fields 



In t 
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ory |32 
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izing the formalism underlying the Andersen-Parrinello-Rahman the- 



] we shall focus on the situation where a system of particles is 
contained in a periodic simulation cell. The formalism we describe can easily be adapted 
to study the deformation of a material confined between walls: to do so, it is sufficient to 
embed the whole system-confined particles plus walls- in the cell and mandate that the 
particles constituting the walls affinely follow its deformation. 

The shape of the simulation cell (which is, by construction, a parallelepiped) is represented 
by the set of d = 2 or 3 Bravais vectors: h = (a,b) or h = (a,b,c); its volume is V = 
det(h) = \h\. We consider a system of N particles, with positions r = {rA- in real space. 
The interaction potential U(r,h) depends on the positions of the particles but also on the 
shape of the simulation cell which enforces boundary conditions. 

"Macroscopic" deformations of the sample are performed by changing the Bravais vec- 
tors. Since we are concerned with variations of the local minimum around some reference 
configuration, we will often use a reference confi gur ation h of the cell and compare it with 
a current value h. Following Ray and Rahman, ^3, [3] we introduce a transformation 
of particle coordinates which maps any vector r onto a cubic reference cell: 

L — h.s , 

with s a G [—0.5,0.5]. If we change the cell coordinate from h to h and require that all 
particles affinely follow the deformation of the cell, any particle at point r is mapped onto r = 
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h.h .r. We denote F = h.h which, in the usual language of elasticity, is the deformation 

== - nn == ■ 

gradient tensor. Once the reference frame h is specified, any configuration ({rj},/i) 

of the system can be parameterized by a pair ({tj},F), with the convention that: r» = 
F.fj. In this parameterization, changes in F_ correspond to affine transformations of all the 
particles following the cell shape, while changes in r$ correspond to the non-affine part of 
the displacement of the particles. 

At zero temperature, an infinitesimal deformation of the system is often performed in 
two steps. First, starting from a local minimum {r { } at h, the particle coordinates affinely 
follow the change of the cell coordinate from h to h. The {f^} remain constant. The 
real-space position of particle i is thus mapped from r_ { onto F^.f^. Second, the particles 
are allowed to relax to the nearest equilibrium position, h being fixed. They reach new 
positions {r^} which differ in general from {F_.ti\- The non-affine part of the deformation 
is then characterized by the displacements as viewed in the reference frame, {-F - .r,}. For 
small displacements, the particles continuously follow changes of local minima. The real 
space positions of the particles at equilibrium are thus a continuous function of h (on some 
interval of strains), and we could denote these equilibria as {r^h)}. Likewise, the continuous 
changes of local minima are most readily studied by monitoring changes in the reference co- 
ordinates: {tiitk) — Er 1 -Liitk)}- 



1.3 Equilibrium trajectories 

By definition, elastic constants are second order derivatives of the energy with respect to 
strain. Since strain is characterized by second order tensors, the elastic constants are fourth 
order tensors. Before presenting the details of the tensorial formalism, however, it seemed 
pedagogically sounder to us to first consider the simple situation where the shape of the cell 
can be parameterized by a single degree of freedom. 

Suppose then that we prescribe the tensor hi^f) as a function of a scalar parameter 7. For 
varying 7, so long as the local minimum does not vanish, the system follows a continuous 
trajectory in configuration space as illustrated on figure ^ Given a reference cell h = h^o) 
at 70, the energy functional can be written either as a function of r and 7: U(r, 7); or as 
a function of r and 7: U(r, 7) = U(F_(j).r, 7). We introduce the notation U to emphasize 
that-contrarily to W-this function is defined after a choice of reference cell with Bravais 
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matrix h. Changing 7 for fixed {f 4 } corresponds to performing an affine strain of the whole 
system-the particles and the boundary. 

When particles are constrained to follow deformation-induced changes of a local min- 
imum, their real-space positions {^(7)} and the corresponding non-affine displacements 
{7^(7)} are now one-parameter functions of 7. An equation of motion for the non-affine 
displacement fields derives from a straightforward application of the implicit function the- 
orem: The trajectory is specified by the condition that the system is always at mechanical 
equilibrium: 



/, = - 



dU 



or, — 







(1) 



Note that the second equality is a property of the point derivatives of any observable of the 
form, ^({7^ = F.fj}): 

(2) 



OA _dA 



and is derived in appendix A. Differentiating the condition Jjr 



leads to: 



d 2 U Vr< 



d 2 U 







once with respect to 7 



(3) 



where the symbol D is introduced to indicate derivatives which are taken under the con- 
straints of mechanical equilibrium. This equation is formally valid for all 7 7^ 70, even 
though it will primarily be used in the limit 7 — > 70 (i.e. h—*K). 

70, equation (J3J) involves the Hessian: 



Using equation (J2J), we see that in the limit 7 

d 2 U 



I(7o) = 

and the field of virtual forces 



7^70 



d 2 U 



(4) 



7^70 



l(7o) 



d 2 u 



(5) 



7^70 



In semi-condensed notation, equation © reads: 



= P7 



(6) 



7^70 



In order to solve equation © we need to take care to eliminate the zero modes of the 
Hessian. Generically there are d zero modes for a system in d dimension: those correspond 
to translation invariance; their existence indicates that a solution to equation © can only 
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FIG. 2: The fields Hj a (left) and T>fi a /T>j (right) as defined in the text for a typical atomistic 
system. The deformation mode is simple shear. Note the random character of Hj Q and the strong 
correlations in Vri a /T)^ 

be defined up to global translations of the particles. There are no zero modes associated 
with rotations because the geometry of the simulation cell (a parallelepiped) breaks the 
invariance of the problem under global rotations of the particles (at fixed cell boundaries). 
The zero modes can thus be eliminated in the standard way by subtracting off the projection 
onto uniform translations from any particular solution. Up to this invariance, the solution 
of (EI) reads: 

= H" 1 .S . (7) 



Vr 



7^70 



The numerical practice which consist in, first affmely deforming the simulation cell, then 
letting the system relax toward the displaced minimum, provides a direct physical interpre- 
tation for this equation. 3 is the field of forces which would result from an elementary affine 
deformation of all the particles. Equation (J2J) shows that the non-afnne displacements are 
just the linear response of the system to these extra forces. 

Examples for the fields S and Dr/Dj in simple shear and pure compression are shown 
figure El and El (The details of the numerical simulation will be given in section 1.5.) We 
observe that the short-range randomness of the vector S contrasts with the large vortex- 
like structures displayed by the non-afnne "velocity" fields Di_jT>^. We will later use the 
apparent disorder of the field S to construct a statistical treatment of the contributions of 
these non-afnne displacement to the elastic constants. For now, let us move on to study how 
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FIG. 3: The fields Hj Q (left) and T>ri a /T>^ (right) in pure compression. 



these non-affme fields enter microscopic equations for the elastic constants. 

Upon deformation, the first and second derivative of the energy Uiff) = ^(£(7), 7) are 
related to the components of the stress and the elastic constant respectively. We will later 
enter a full tensorial derivation of the equations to specify this relation. For now, we can 
simply differentiate U{^f) with respect to 7. At first order, we obtain a projection of the 
stress: 

Til J r)7°J T)r Rl°J Rl°J 

(8) 



du 

$7 



du 



D7 dr P7 

where the second equality holds because of mechanical equilibrium. The partial derivative 
which appears at the rightmost side of these equalities corresponds to the projection of 
stress which would be inferred from the assumption that all the particles undergo affine 
displacements. Since the total derivative of the energy as a function of 7 (while enforcing 
the condition of mechanical equilibrium) is identical to the partial derivative, it means that 
non-affine displacements do not contribute to the first derivative of energy with respect to 
7, i.e. the projected stress. 

The situation, however, is different when the second derivative of the energy is considered. 
Indeed, we have: 



V 2 U 



d 2 U d 2 U Vr 

+ 



d 2 U „ 



£>7 2 d'-f 2 djdr P7 d'-f 2 — 



S.H '.S 



(9) 



The first term in the right hand side of this equation, corresponds to energy changes 
associated with strictly affine displacements of the particles. It is the Born approximation for 
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the second derivative of the energy with respect to 7. The last term introduces a correction to 
the Born approximation which results from the existence of non-affine displacements. Since 
H 1 is positive definite, the corrective term is negative: non-affine displacements reduce 
the second order derivative of the energy from its Born approximation. This is another 
statement of the idea that non-affine displacement are associated with further minimization 
of the energy of the system from a tentative homogeneous deformation. This property holds 
for any mode of deformation h('j). It means, for example, that the non-affine corrections to 
the shear and compression moduli must be negative. 



1.4 Tensorial forms 



In the preceding section, we considered the non-affine displacements in response to some 
prescribed mode of deformation parameterized by a scalar parameter, 7. We now proceed 
to study the tensorial form of the above equations: derivatives of the energy functional with 
respect to the components of the strain tensor will provide analytical expressions for stress 
and elastic constants. 

The calculations can be simplified by noting that the number of independent components 
of the strain tensor is reduced by symmetry. In particular, a useful symmetry property 
holds under the general assumption that 
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42j the interaction potential can 



be written as a function of the full set of distances between pairs of interacting particles: 



W(r,/0 = U ({? 



,-|}), where the index ij runs overs all pairs of particles, and for each 



pair, r_ii = r_j ~ Li modulo the periodic boundary conditions. (Note that for a periodic 
simulation cell, the distance between two points depends on their positions relative to the 
boundaries of the cell: the way distances are calculated thus depends on the shape of the 
cell, hence, on h.) We stress that the ability to write the energy functional in this way 
is not restricted to potentials involving only pair interactions: it also holds for e.g. three 
body interactions which depends on bond angles, as used for silicon or embedded 
atom potentials 0, as used for metals. This formalism does not apply, however, to 
situations when rotational degrees of freedom must be taken into account as, for instance, 
in Cosserat apaches t ^ — HHQQorto .ode! atopic object, 



li k e ne mat ics. Q We expeet a si.ni.ar toal t—t to be possibie in tbe 8 e situations, 
but it will require using a slightly more general formalism that we do not wish to address 
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here. 

Let us consider a reference configuration h and a current configuration h. Suppose that 
r and r are the difference between the positions of two particles in both these systems of 
coordinates. They are related by r — F_ . r, whence: 

r 2 = 2 r T .ri.r , (10) 



2 r .rj.r 



with the Green-Saint Venant strain tensor: |35 

1 



36 



37] 



n 



(ii) 



We want to write the energy functional after the change of variable where any config- 
uration of the system is mapped onto ({f^} = F_ -Ln^O- Under the general 
assumption that the potential energy can be written as a function of the distances {r^} 
only, equation ([TO]) can be used to write: 

U(r,h) 



U 



ri + 2 rJj-V-tij 



-ij 



Where it appears that the energy of the system is a function of {tij} an d r\ only. [l| 



41 



40 



42| The notation U is introduced here to emphasize that this functional of {tij} an d f] is 



defined after a choice of reference cell with Bravais matrix h (whereas the functional U(r,h) 
does not depend on this choice). A choice of h being made, however, the energy functional 
depends on the current cell coordinates only ma the symmetric tensor r], but not loathe whole 
tensor h. We recover the separation of coordinates: changing 77 for fixed {r^} corresponds to 
performing an affine strain of the whole system-the particles and the boundary; changing 
{r^} in the reference cell corresponds to performing non- affine displacements of the particles. 

We are now in a position to derive equations of motion for non-affine displacement fields 
in a fully tensorial form. In any configuration, the force acting on particle i is: 



i=-^L(fe}.fi>=-||,(fe}.2)-r 



and mechanical equilibrium reads: 

y . dU 



dU 



dh It; 



(12) 



(13) 



The derivatives are taken here for fixed cell coordinates. The equation of motion for {fj}, 
is now obtained by differentiation of (jl3|) with respect to the components of rj\ 



d 2 u 



+ 



d 2 U 







(14) 
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with the usual summation convention of repeated (Greek and Latin) indices. As above, 
the symbol T> denotes partial derivatives with respect to some tensorial components while 
enforcing mechanical equilibrium ([13)1 . In the limit r\ — > 0, equation (fTlj) involves the 
Hessian: 



and the field of third order tensors, 



h=h 



d 2 u 

dr 4 dr 4 



(15) 



T)->0 



d 2 u 



dr_t drj l 



(16) 



In semi-condensed notation, equation (fT^j) reads: 

Vr 



(17) 



Provided standard caution in the inversion of the Hessian, we then get, up to translation 
modes: 

H- 1 ^ . (18) 



which is the non-affine displacement field in response to elementary deformation along rj KX . 

To get further insight into the correction term in equation ()24|). and on the physical 
interpretation of 3, let us write: 

<9f. 



dr), 



(19) 



77— *o 



This equality comes after differentiation of equation (|12j) with respect to 77 and use of me- 
chanical equilibrium (equation (fTT^) ). As in our pedagogical example of one-parameter strain, 
the field S KX , can be interpreted as the forces which would result from an elementary affine 
displacement of all particles in the strain direction r) KX . (More specifically, At7 kx H kx , with 
no summation on the greek indices, is the force resulting from a small affine transformation 
by Ar] KX ; S kx is the tangent direction to the changes of forces upon affine transformations.) 
We see that Sj AX can be seen as the fluctuation of the force /. in response to an elementary 
strain. This interpretation emphasizes the random character of field S KX : local forces are, 
by definition zero at mechanical equilibrium; however, their variation under an elementary 
strain depends on the configuration of the particles with which particle i interacts; in partic- 
ular, KX should be zero if the conformation of particles surrounding i is symmetric-which 
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is the case of a Bravais crystal. H KX is thus a measure of the local asymmetry of particle 
configurations. 

Equation (|18|) further states that the non-afline displacement is nothing but the linear 
response of the particles to these extra forces. 

Let us now derive microscopic equations for the stress and elastic constants. Because 
the total potential energy is a function of 77 and r only it is sufficient to characterize the 
elastic response of a material using the derivatives of the potential energy with respect to 
the components of 77. We note, however, that the usual definition of the Cauchy stress and 
elastic stiffnesses involves derivatives with respect to the components of the deformation 
gradient tensor F_. Derivatives with respect to the components of r\ provide the so-called 
thermodynamic tension and elastic constants (or also thermodynamic stiffnesses). The def- 
initions of thermodynamic tension, Cauchy stress, elastic constants, elastic stiffnesses and 
their relation are given in appendix A. 

By definition, the thermodynamic tension is the derivative of the energy functional with 
respect to the components of r\: 

_ 1 VU 

= ~v?>n 

Since 77 is a symmetric tensor, we can conclude without further examination that thermo- 
dynamic tension is also symmetric. Using mechanical equilibrium, we next have: 

I VU 1 dti 
t — — — — = — — — — , 20 

= v T>n v dn 

As in our pedagogical example of one-parameter strain, we see here that the thermodynamic 
tension is equal to the partial derivative of U (here, with respect to 77). This partial derivative 
corresponds to changes in energy during a strictly affine displacement of the particles: we 
see that the existence of non-affine displacement fields does not appear in the expression for 
the thermodynamic tension. Expression f|20[) is valid for any finite strain: this will allow 
us to later differentiate this expression a second time and obtain elastic constants. Often, 
though, we will use it in the limit 77 — > 0. 

Note that this expression for the thermodynamic tension provides another interpretation 
of the field S. We can write: 



(21) 

77=0 
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from the definitions of t and S. This expression is to compare with equation ()19|). We see 
that Hj KX can be seen as a fluctuation of stress (or tension ) in response to an elementary 
displacement of particle i. This interpretation was noted in 14j. 

The elastic constants are second derivatives of the energy functional with respect to the 
Green-Saint Venant strain tensor: 

1 V 2 U 



■0KX 



(22) 



The derivative is here "total" since this is the second derivative of the energy following 
deformation-induce changes of a minimum. As we can permute the order of derivatives, 



elastic constants verify C a p KX = C Kxa /3, and since r\ is symmetric, C a p KX = C3, 



CtKX 



To obtain an analytical expression for the elastic constants, it suffices to differentiate ex- 
pression ()2())1 once and take the limit 77 — >• afterwards: 



d 2 U d 2 U Vr, 

+ 



Using equation (fTSj) and the definition of 3 a/3 , it then comes: 



(23) 



??->o 



a 



d 2 U 



afiKX 



V \ dr} a pdr] KX 



" I I 



(24) 



We recognize the first term in equation ()24)1 as the Born approximation C^. The con- 
traction of the inverse of the Hessian on components of S a/3 provides the correction terms. 



1.5 Special case: Isotropic material in 2D 

Here, we illustrate these ideas with numerical simulations of a two-dimensional bidisperse 
mixture of particles interacting through a shifted Lennard- Jones potential. Particle sizes 
r$ = r^sin ^5 /sin f and a number ratio N^/Ns = 1+ 4 are used to prevent crystallization. 
The system is prepared via an initial quench from an infinite temperature state. Further de- 
tails regarding the numerical protocols will be found in our forthcoming dedicated numerical 
study |2q . 

The interaction potential is pair- wise additive: W({r i:) }) = (r^). In this case, 

formal expressions for the fields KX and the Born approximation for the elastic constants 
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can be related directly to the derivatives of Vy. These expressions ((J7TJ) and (|72|l) are derived 
in appendix B. We just recall here equation (J72*Jl : 

^a/3Kx = ^ ^ Ci 3 ~ tij) r ij n ij n ij n ij n }j (25) 

ij 



where we have introduced the normalized vector between pairs of particles: n i5 - = — , as 



r ■ 

-t- 

well as the bond tensions and stiffnesses: 



tin = — — - and c. 



13 drij 13 <)rf j 

The material is expected to be isotropic so that the elastic constants take only two 
independent values under permutations of indices: C a p KX = A 5 a p 5 KX + /i (5 ax 6p K + 5 aK 5p x ), 
which define the Lame constants, A and /x. Using equation f)25j) we observe that the Born 
term is identical for the two Lame constants A and \i. Introducing the notation = 
(cos%, sin%), we find: 



A Born = ^Born = I ^ ( r< . Q . _ tij ) Tij C OS 2 (%) sin 2 ( 



To calculate the Lame constants, it suffices to consider two modes of deformation and 
the associated non-affme fields. We thus consider S s = = S yx and 3^ = (S xx + 'E yy )/2 
which are associated with pure shear and pure compression respectively. Plots of S s and 
H c , and their corresponding displacement fields were given on figure 121 and El respectively. 
Pure shear grants access to the sum: C xyxy + C xyyx + C yxyx + C yxxy — 4 /i. Pure compression 
grants access to the sum: C xxxx + C xxyy + C yyxx + C yyyy = 4 K, where K = A + \i is-in two 
dimensions- the compression modulus. With the fields H s and S c , we can directly construct 
the non-affine correction to /i: \x — — S S .H _1 .S S ; and to K: K = — H C .H _1 .H C . We recall 
that since the Hessian is positive definite, these corrections are necessarily negative. We thus 
have in all generality that for an isotropic material, \i < /i Born and K < K BoTn . Using these 
fields, we find in this sample, for the shear and compression moduli: /i Born ~ 125, jl ~ 86, 
whence \x ~ 39 and K Born ~ 250, K ~ 14, whence K ~ 236. All moduli are reported in 
dimensionless stress units. These values compare to the highest density systems studied by 



Tanguy et al 



10 



111 



We note that the correction term to the compression modulus appears to be small. This 
can be understood to arise from a simple property. Let us consider the case when the po- 
tential is pairwise additive, and is homogeneous in the sense that the force between each 
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pair of particles is a homogeneous function of their distance: F(nrij) = K a i^(r»j). A pure 
compression corresponds to a global scaling of all distances: if the forces are homogeneous, 
a global scaling of the distances preserves a state of mechanical equilibrium. In other word, 
whenever the forces are homogeneous functions of the distances, the fluctuating force field 
associated with pure compression, S c , vanishes identically. Hence, the correction to the Born 
approximation vanishes for the compression modulus: K = A + \i — K Born = A Born + /i Born . 
For a compressed Lennard- Jones system, if the total energy is dominated by the pairs of 
particles which are at close distance, the interaction is thus dominated by the repulsive 
power-law divergence of the potential. Since this power-law leads to forces which are homo- 
geneous functions of the distances, this repulsive part of the Lennard- Jones potential leads 
to vanishing corrections to the Born approximation. As a result the overall amplitude of the 
field S c is small, hence the small correction to K. In particular, this is consistent with the 
data from Tanguy et al jloL Q] and explains why the corrections induce important changes 
in both A and [i but smaller changes in their sum. 



2. NORMAL MODE DECOMPOSITION 



In expression (|24|). the correction to elastic constants involves a contraction of the Hessian 
on the fields SLg. Since the Hessian is positive definite- we are perturbing around a local 
minimum-it acts as a scalar product in expressions such as: E aj g.H -1 .S . This form of the 
corrections to elasticity suggests that further insight may be gained from a normal mode 
decomposition of the fields S KX . 

We denote the eigenvectors of the Hessian by and the associated eigenvalues by X p . 
We introduce a particle mass to correctly scale eigenvalues and eigenfrequencies and write 
X p = mujp. (In our numerical simulations, the mass is taken to be unity.) The vector S KX 
is written as: 

s =v s * p 



with S PiKX = H K ^.\I[_ P . Using this decomposition and equation (JTHJ), the non-affine displace- 

J2=f^* p ; (26) 



ment fields read: 

Vr 
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the non-affine contribution to elasticity reads: 

1=1 it 1 r-i \^ "P,a/3 "PAX 



X„ 



As we saw from figures 121 and El the fields S display a very random character. The 
randomness of the fields 3 can be further assessed by studying fluctuations of the weights of 
their normal mode decomposition. To examine this question, we report, in figure EJ all values 
of H P)W , with u = s, c corresponding to pure shear deformation and pure compression. These 
quantities are obtained by calculating the fields H u and performing a full diagonalization of 
the Hessian. We see that within a frequency shell tu p G [lu,lu + da;], the values of are 
scattered, so that they can indeed be interpreted as random variables. On the scatter plot, 
it is apparent that the distributions of S p u are symmetric with respect to the horizontal 
axis, hence that the H pu 's have zero mean: this property is enforced by symmetry since 
each of the eigenvectors of the Hessian is defined up to a sign convention. The width of the 
distributions of seem to depend smoothly on uo and increase with oj for the most part. 
At high frequencies the scatter plots thin out because there are fewer and fewer eigenvectors 
to sample the density of states close to its upper cut-off. 

On the basis of this observation, we make the assumption that the random fields, S KX , 
and in particular their projections on normal modes, S P)WX , are self-averaging quantities. We 
thus introduce the correlators on frequency shells: 

Ta/JrextV) = ( s p,a/8 s p,«x)w p e[w,t«H-du;] (28) 

In this definition, the average is performed for all the projections of 3 Q/3 and 3 KX on eigen- 
vectors with eigenfrequency lu p G [lu,lu + da;]. The assumption that S PjKX are self- averaging 
means that the quantities T a p KX (uj) are expected to converge toward well-defined functions 
of uo in the thermodynamic limit: either when a large ensemble of systems of a given size is 
considered, or when one single system of a large size is considered. 

From their definition, we observe that functions T a p KX (uj) verify the same symmetries as 
the elastic constants, namely: 

r ^ aflnx ^ KX a P ^ P aK X (29) 

In the thermodynamic limit, we can finally rewrite equation (|24jl as: 

r ^yBorn f°° j, , P( U ) Fg(3 KX (uj) 

^apK X - O q/3kx - / do; —— 2 , (dUj 

In IIIUJ 
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FIG. 4: A scatter plot of the quantities H Pjtt versus the corresponding eigenfrequencies of the 
Hessian (with u = s in shear (top) and u = c in compression (bottom)). B. pu is the projection of 
the field S u onto the eigenmode with frequency u p . This set of data has been obtained using 
one typical system of size L = 50. Each point makes a contribution to the non-affine corrections 
to the elastic constants as described in the text. 

where p is the density of states This equation is a sum rule: it relates elastic constants 
to the integral of a correlator between microscopic fields. 

Like the elastic constants, T a p KX are tensorial quantities. Therefore, if the material is 
isotropic-given the symmetry property (J2HJ)- these functions can take only two independent 
values corresponding to the two Lame constants A and /i: 

r Q/ 3 KX (w) = T A (u;) 5 a p S KX + T^(cu) (5 ax 5p K + 5 aK 5p x ) 

In this case, the two correlators, T\ and provide a complete description of the non-affine 
corrections to the Lame constants according to equation (jHOj) . 

To illustrate this discussion, we return to our two-dimensional numerical study of two 
modes of deformation, simple shear and pure compression. This protocol suffices to measure 
the two Lame constants of our system. The vector fields S s and S c grant direct access to 
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the quantities T u = (2^ u ) for u — s,c. Elementary algebra, similar to the calculations of 
section 1.5, shows that these quantities also verify: = T s and IV = T x + T M = T c . 

We report, in figureEl the density of states p(co), the correlators T s , T c , and the quantities 
C u (u) = p{oS) T u (oj) / (mu) 2 ) (with u — s,c) which measure the total contribution of the shell 
[u,u + du} to the non-affine corrections. Three system sizes have been used in this plot: 
L = 20, 25, and 30, where L is the length of the square simulation cell in dimensionless units. 
For each size, an ensemble of 20 system was used to provide statistical accuracy. 

The functions T s = and T c = Tk measure the variance of the distributions of 2 PjS 
and 2 PjC respectively, which appear on the scatter plot of figure 13] We see that that these 
functions are roughly increasing as was suggested from the observation of the scatter plot. 
Important fluctuations of V ^ and IV at high frequencies result from the lack of statistical 
representation in the region where the density of states vanishes and where the scatter 
plots thin out. For their main part, however, the curves lie on top of one another for the 
three system sizes we have used: this is evidence that the continuous limit has already 
heen reached. This obs e r va tlo n is con Slstent with t he _ o f Ta„ guy et a, Q Q 
since the values of the elastic constants they measured did not vary much for the different 
sizes they studied, from L ~ 20 to L = 150. We believe that the good convergence of 
these measurements results from the decorrelation of the fluctuating forces 2^ at very short 
distances: because the correlation length is short, the finite sizes of the systems considered 
here provide enough sampling of the distribution of 2^ 

The density of states, and the T's exhibit rather complicated functional forms. In contrast, 
it is striking to us that for simple shear and pure compression C u (u) seems rather featureless; 
particularly in comparison with both the density of states and the functions T u (lj). We have 
observed this for different numerical models and different densities. The monotonous 
decay of C m (cj) and Ck{w) seems to indicate that a simple physical mechanism-like a trans- 
fer of elasticity from low to high frequency phonons-governs the non-affine corrections. This 
reminds us of an analysis by Radjai and Roux of quasi-static deformation of granular mate- 
rials ^] in which the authors suggested that the non-affine velocity field could exhibit statis- 
tical properties similar to turbulence. This sound like an appealing suggestion, although, at 
present, we have no evidence to further test the existence of turbulent-like energy transfers 
from large to small scales. 
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FIG. 5: Ensemble averaged contribution to non-affine elasticity from individual bands, [uj,uj + duj], 
as described in the text. Results were obtained as ensemble averages over 20 systems of sizes L 2 
with L: 20 - circles (black), 25 - squares (red), and 30 - diamonds (green). Convergence of the 
curves demonstrates that finite size effects are small, (a) Density of states. (b,c) Magnitude of 
"noise" field, {SS)(u), for simple shear (b) and compression (c). (d,e) Net contribution to elasticity, 



, of each frequency band for shear (d) and compression (e) . 
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3. VISCO-ELASTICITY 



In the above formalism, we assumed that deformation was very slowly forced upon a piece 
of material, so that the relaxation of the system toward an inherent structure was always 
faster than the changes induced by the external drive. This guaranteed that the system 
followed an energy minimum in configuration space. Our analysis, however, is essentially 
based on a second order expansion of the energy functional: it assumes small displacements, 
but may, in principle, apply to situations when the timescale of the external drive becomes 
relevant. 

We consider here oscillatory perturbations of small amplitude, and assume that the sys- 
tem remains in the linear response regime around some energy minimum, far away from 
any catastrophic event. We expect this situation to be of relevance to questions such as e.g. 



acoustic damping in granular materials [17|, |18|, |19|, |20| or viscoelasticity in dense emulsions 
and foams Near a minimum, but not exactly at mechanical equilibrium, the 

particles are subjected to non-zero forces: an oscillatory external forcing at finite rate thus 
works on the system of particles. Energy is injected into the system and is lost in dissipative 
mechanisms at the microscopic scale (friction between grains or viscous dissipation between 
bubbles). 

We will follow the usual practice of characterizing energy transfers in the linear response 
regime using the components of the complex stress response, that is the storage and loss 
moduli, G' and G". j^jj We insist that we focus here on features of G' and G" that arise 
merely from the existence of non-affine displacement fields around a single minimum. In 
an experiment-think, for example, of a foam-it is unclear whether sufficiently small strain 
amplitudes can be achieved so as to stay within the basin of attraction of a single equilibrium 
configuration. Thus measured values of G' and G" would also receive contributions from 
plasticity in addition to the microscopic visco-elasticity studied here: transitions between 
distinct inherent structures would become a relevant mechanism of energy dissipation [3, 
0]. It's only after a careful study of the contribution of linear response that we will be in 
a position to characterize the relative importance of these various dissipative mechanisms. 
The present work is a preliminary step in this direction. 
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3.1 Damped equations of motion 



In order to work at finite strain rates, we need to use Newton's equations coupled with 
the deformation of our simulation cell. In order to prevent the system from heating up, a 
microscopic mechanism must be specified for energy dissipation. Dissipative mechanisms 
vary from system to system. In granular packings, dissipation involves the deformation 
of asperities at contact between grains. In foams, dissipation involves fluid flow inside a 
liquid phase: it is modeled by Durian 0, |2| by viscous terms involving velocity differences 
between neighboring bubbles. In numerical simulations of structural glasses and supercooled 
liquids, the coupling with an external bath is modeled either by a Nose-Hoover 



thermostat, by use of a Lagrange multiplier (59, 



60j], or by a Berendsen thermostat [61 




It thus appears, after examination of these different systems, that energy dissipation often 
arises as some sort of viscous damping at the microscopic level. We will thus limit our 
present discussion to the situation where dissipation enters via a viscous drag force applied 
on individual particles. Below, we give some rule-of-thumb estimates for what the value of 
this viscosity might be in the case of dense emulsions. 

In the SLLOD algorithm, the viscous term is generally taken to be proportional to the 
peculiar velocity, 6JJ that is the velocity r_ t of the particle minus its velocity corresponding 
to the affine flow ^(rj = F.f^ = F^-F^ 1 .^. Similarly, in foams or granular materials, the 
viscous damping does not depend on the velocity in real space but on the difference between 
the particles' velocity and their surrounding environment. We thus first write the equations 
of motion as: 

m li = Li~ v (n-n(ju)) ( 31 ) 

In term of {ti\, it follows that: 

df„- 



d 2 f- ■■ • dr ■ 

mF.—^- = f -mF.r i -2mF.^-uF.^ 
= dt 2 =~ l = dt = dt 



(32) 



We see here that a term mF.fj arises which depends on the position of the particle in 
the simulation cell. Another term of similar flavor would be v F_.r_{. it does not appear 
in equation ()32|) because it has already been eliminated by the above assumption that 
the viscous dissipation applies to the peculiar velocities. Such terms introduce a spatially 
dependent coupling between particular motion and the cell coordinates: they cannot be 
allowed to enter the equations of motion if we want to enforce translation invariance. This 
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was originally noted by Andersen 



32| and later by Ray and Rahman |36j. Within the 



3- 



formalism of Souza and Martins [63[, an equation of motion with no spatially dependent 

terms can be derived systematically from a Lagrangian. It relies on the idea that the distance 

traveled by a particle should be calculated while removing a component of rigid rotation. 

We do not wish to detail this formalism here, but simply recall that adding appropriate 

terms in the Lagrangian of the system (particles plus simulation cell) permits eliminating 

space-dependent terms from equation (|32j). Following these authors, we eliminate the term 

mF_.r 4 and focus on translation invariant equations of motion for an athermal system in 

SLLOD dynamics. These equations read: 

d 2 r- • dr ■ dr- 

= dt 2 - l = dt = dt ( J 

We want to account for the dynamics of the system when it is submitted to small amounts 
of strain \\F_ — 1\\ ~ e << 1. To do so, we write a perturbative expansion of equation ()33|) 
around a known equilibrium state {r^}. This expansion is written in terms of the displace- 
ments {Xj(t) = fj(t) — fj}, which are also of order e. At first order in e, the perturbed 
equations of motion read: 

d 2 Xt x + . v <^L 

dt 2 dr_A ~ J dt] = dt 



= -J=±.xj + ^ =1 -V-v^ ± (34) 



Using the definition of H and 3 (equations (fT3j) and (fT6*j) ) it follows that: 

m ^ + ^ + ^'^ = ^^ (35) 
This equation governs the linear response of our system around a single minimum. Starting 
from this equation, various limits allow us to recover more usual expressions. For example, 
taking v = and canceling the term on the right hand side yields the equation governing 
the vibration modes of the system; canceling m and v yields the equations which governs 
quasi-static response and defines the non-affine displacement fields. 

To solve equation (|33j) . we perform two transformations: a Fourier transform of the 
time domain, and a normal mode decomposition of the displacement field. The Fourier 
components 5?j are the responses to perturbations of the form r](t) = rjsin(ut). They are 
further decomposed onto normal modes as: x_i(uj) = x p (u)^ and x_ p (uj) = Xi(u)) With 
these notations, we find: 

— m uj 2 x p + i uj ux p + muj 2 Xp = H pAX rj KX (36) 
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which is solved by: 

~ . "PAX Vk-x 



muj 2 — muj 2 — iujv 



v 

The perturbed stress is then obtained by writing at first order in e: 



1 / d 2 U d 2 U 



V \dr) a pdr] KX drj a pdti' l y 

Using (j37j), the complex stress response in the frequency domain thus reads: 

p 

= + i v — f^P^x — ~ {u) 

a/JKX l*X\ J y muJ 2 _ muJ 2_ luJhl ' K X\ I 

p P 

It can thus be recast in the form: 

At a p(uj) = G al3KX (uj)rj KX (uj) 

with: 



v 



(37) 



n (, ,\ r^Born | 1 "P,a/3 "PAX /qo\ 



Taking the thermodynamic limit, we can furthermore introduce the functions r Q ,^ KX (u;) 
defined in equation (J28)) . Equation (|38|) can then be written as an integral in frequency 
domain: 



<W") = CXl + I <H - T^l^T... • (39) 



OO 



p(^p) r 

q/3kx (^p) 



1/3kx ' p muj 2 — muj 2 — iu) v 



It is then an easy task to extract the real and imaginary part of the stress response. Observe 
that we recover the Born approximation and the true elastic constants in the high and low 
frequency limits respectively. 



3.2 Overdamped systems and the relaxation spectrum 



The response moduli are often studied in the 
remains small compared to the elastic energy. 



overdamped' limit, when the kinetic energy 



5 it ] In this case, the term muo 2 is negligible 
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before mu p and equation can be given a simpler form. To do so, we introduce the 
timescales: 



and write equation (jSUjl as: 



v v 
p mull A„ 



1 Z" 00 .. 1-iuTr, 



<W<") = CXc " o / d In r p — — -f H(r p 



1 Z" 00 icjr„ 
C «^x+o / dhiTp y g(r p ) 

2 y_oo 1 + 2 W T v 



where 



H a p KX (r) — \/ —— p ( W— J r a/3KX ( W— ) (40) 



The complex modulus G a p KX is thus naturally written as a sum of Maxwell elements. By 
definition, J^jJ the function H a p KX is the relaxation spectrum of the system. It is the prod- 
uct of the density of states with the correlator T a p KX . The specific forms of this relation 
associated with a given experimental geometry -pure shear or pure compression-derive triv- 
ially from these tensorial expressions. We plot the shear relaxation spectrum for our model 
system in figure again for each of the three ensembles of systems with lengths 20,25, and 
30, demonstrating the system size independence of the results. 

Since the true microscopic mechanisms of dissipation in real structural glasses remain 
somewhat poorly understood, we focus here on the case of a dense emulsion or foam, where 
the system is overdamped by construction. To make a semi-quantitative comparison with 
experiments on such systems, we need to give a rough estimate for the values of the coefficient 
of viscosity, v, and the overall stiffness scale, A, in our model. A truly quantitative connection 
with experiment is somewhat beyond the scope of the present work and will be explored in 
our dedicated numerical study [25j . We should also remind the reader that the data presented 
in this work for illustrative purposes is taken from a compressed system with Lennard- Jones 
interactions... clearly not sufficient for a fully quantitative comparison. However, the hope 
is that there is a good deal of qualitative similarity in the structure of the Hessian matrix 
for these various classes of disordered systems, and we proceed to dimensionalize A and v in 
order to make a semi-quantitative comparison with experiments on packings of emulsified 
droplets. 

A series of experiments, simulations, and theoretical works on such a system was carried 



out by Mason and coworkers 
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23j. The surface tension of the droplets was found to 
26 



100 




FIG. 6: The relaxation spectrum (as denned in the text) for simple shear. As above, each of the 
three curves corresponds to a 20 fold ensemble of square systems of length, L: 20 - circles (black), 
25 - squares (red), and 30 - diamonds (green). 

be about a ~ 10 dyne/cm, which, for well-compressed systems, should set the average scale 
of the interparticle stiffnesses. For crystalline systems, the scale of interparticle stiffness 
is characterized by the maximum eigenvalue of the system (corresponding to eigenvectors 
which are plane waves with periodicity equal to one reciprocal lattice vector), thus we scale 
the eigenvalues of the Hessian matrix of our system such that the maximum eigenvalues 
are roughly a. Those maximum eigenvalues of the Hessian matrix, as can be read off from 
the density of states, are about 5000 in our dimensionless units of stiffness, so we must 
have A = pr = = 2.10~ 3 dyne/cm. This sets the scales of the second derivatives of 
the potential. The order of magnitude of elastic constants then comes out by dividing A 
by the average interparticule distance, that is by the diameter of droplets in the case of 
an emulsion. Taking a droplet radius of a micron, this sets the scale of elastic constants 
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as [C] = 10 dyne/cm 2 . Our measurement of elastic constants in the range of 100, thus 



corresponds to values o 



order 10 3 dyne /cm 2 



Following Durian |55j], we ought to choose a value for v such that the drag force, F^ Tag — 
u8v, between two droplets sliding past each other with relative velocity Sv generates a stress 
in the film which is equal to the strain rate in the film times the effective viscosity, so we 
must have: v = Fd rag /6v = VeS^r where R is the typical lateral extent of the film (which, 
following Durian, we take to be of the order of the droplet radius itself), and / is the width 
of the film gap. We can now dimensionalize our units of time, [T] = u/X, using the estimate 
of Liu et. at. for the effective viscosity in the film at the droplet interface of lcp, a droplet 
radius of a micron, and a film thickness of a nanometer, yielding: [T] = 5s. From the plot 
of the relaxation spectrum in figure El we expect the viscoelastic effects and the crossover to 
affine deformation to appear in the range of frequencies from around 1Hz up to about 1kHz 
in these systems. 



4. CONCLUSION 



In this article, we have studied molecular displacements associated with deformation- 
induced, continuous changes of a minimum in the potential energy landscape of an amor- 
phous solid. For small amounts of deformation, the trajectory followed by the particles is 
smooth, and the tangent displacements-analogous to a velocity, with strain playing the role 
of an effective time-involve affine and non-affine fields. Both fields enter at first order in the 
equation of motion of the particles during quasi-static deformation. 

The non-affine fields can be calculated by inverting the Hessian on a fluctuating force field 
S, which is formally a third order tensor. For every component r] a p of the deformation, SLg 
is the field of forces resulting on every particle from the affine motion of its neighborhood. 
The non-affine displacement fields and the corrections they induce on elastic constants, 
can be expressed solely in terms of the Hessian and of various tensorial components of the 
field ELg. We next observe that the field is weakly correlated in space, so that it can 
essentially be considered as a random vector field. This randomness is confirmed by the 
scattered values for the projections of S a/3 on the eigenmodes of the Hessian. 

The normal-mode analysis of the field ELg grants access to the contribution of each 
frequency shell to non-affine corrections to elasticity. Our numerical calculations indicate 
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that the correlation function T a p KX (uj) = (^ p>a j3 ^p,kx)w p e[u,u+Sw] converges toward a smooth 
function in the thermodynamic limit: this means that the contribution of each frequency 
shell to non-affine corrections to elasticity is self-averaging. In our view, the existence of a 
well-defined limit for the elastic properties of amorphous structures, as observed by Tanguy 
et al, rest on this self-averaging property. 

We moreover observed that the contribution of frequency shells seems to be a rather 
simple function of the frequency, as opposed to both T a/3K _ x (uj) and the density of states. 
This simple form of corrections to elasticity suggests that an elementary mechanism of 
transfer of energy between frequency shells may be at work and determines corrections to 
elastic constants. This is reminiscent of the observation by Radjai and Roux of turbulent 
like features of non-affine displacement fields, [^] although we have not determined whether 
turbulent- like scalings arise in our non-affine fields. A better understanding of the mechanism 
underlying this transfer through frequency shells should in principle permit construction of 
approximations for the elastic constants of amorphous materials. 

After studying elastic properties of an amorphous solid in response to quasi-static de- 
formation, we considered the case when the deformation rate is finite. This required us 
to introduce both a molecular mass and a damping term to provide equations of motions 
which correspond to a molecular system in contact with a bath at low or zero temperature. 
We have shown that the visco-elastic response of the solid can be written in the form of 
a sum of elementary damped oscillators. The frequency spectrum-that is the distribution 
of timescales of elementary vibration modes- is directly related to the function T a p KX {uj). 
It thus appears, that a broad spectrum of timescales arise at linear order solely from the 
structure of any particular energy minimum in a high dimensional configuration space. The 
system remains solid since it resides in the neighborhood of a local minimum in configu- 
ration space at all times: this is different from the usual idea that visco-elastic response 
is associated with transitions between minima in the energy landscape 0, H]- m a rea ^ 
system, this type of dissipation may occur simultaneously with other mechanisms: energy 
transport via phonons, dissipation resulting from anharmonic terms, or transitions between 
various energy minima {i.e. true plasticity). The estimation of these various contributions 
and their interplay seem to us an exciting direction for future research. It will require the 
construction of models, close to the experiments, where various dissipative mechanisms can 
be estimated quantitatively. 
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Finally, we stress that the corrections to elastic constants involve l/co? factors as in 
equation (J5UJl : the convergence of the integral in equation (j3T?j) thus depends on the low 
frequency behavior of the functions T a p KX and of the density of states. This becomes a 
problem when the system develops low frequency eigenmodes which have a non-zero scalar 
product with the field S: a quick inspection of equation (J27|) indeed shows that such localized 
low frequency phonons would lead to diverging terms in the non-affine corrections to the Born 
approximation. We can identify two situations when this occurs: (i) A localized eigenvector 
with a vanishing frequency is involved whenever the local minimum occupied by the system 
reaches a catastrophe. 0, Q This property enabled us to obtain universal scalings for 
the elastic moduli close to a shear induced catastrophe, [3^ and observe the divergence of 
the non-affine corrections at these points, (ii) Eigenmodes seem to drift toward the low 
frequency part of the spectrum when the system approaches the unjamming point of Liu 



and Nagel 
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67j A divergence of the non-affine corrections to elasticity could thus 
arise around this point. In this case, we expect the divergence of the non-affine corrections 
to elasticity to control the unjamming transition around the J point. 
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APPENDIX A: STRESS AND ELASTIC CONSTANTS 



Stress and tension 
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681 Since 



Stresses are first derivatives of an energy functional with respect to strains, 
there are different definitions of the strain tensors (77 or F_ or u = F_—l or e = |(w-Ht T )), there 
are different ways to take this derivative. Moreover, this derivative can be Lagrangian or 
Eulerian. We thus have various choices, which leads to a number of possible definitions. Two 
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definitions of the stress are most important: the thermodynamic tension is the Lagrangian 
derivative with respect to 77; the true stress is the Eulerian derivative with respect to the 
deformation gradient tensor F_. 

Let us consider some energy functional, parameterized by the cell coordinates: W(h). 
This functional can represent different objects: for example, it can be the energy ZY({fj}, h) 
for fixed {r 4 }: its strain derivatives were denoted as partial derivatives in the previous 
discussion; it can also be the energy U({Zi{lk)},lk), provided constraints which enforce a 
relation {fj(/Y)}. The formalism developed here does not depend on any specific definition 
of this energy functional but only on the existence of some function W(h). The value of W 
in a given reference configuration is denoted W. 

We saw the important role played by the Green-Saint Venant strain tensor rj: it accounts 
for the mapping of distances after an affine transformation (see equation (|10[l). This property 
permits writing the strain-dependence of the energy functional via 77 only, whenever the 
energy is a function of the set of distances {r^-} between the particles. |4l| By definition, the 



thermodynamic tension 



ia 



3 



0] is conjugate to the Green-Saint Venant strain tensor: 

i - IT («) 

- V dr] 



Where V is the volume of the reference cell. The thermodynamic tension t is defined after 
a choice of a reference configuration: it is a Lagrangian derivative. It identifies, in the 
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continuum limit, with the second Piola-Kirchhoff stress tensor. |3£ 
known, in general, not to have a simple mechanical interpretation: it cannot be interpreted 
as the tensor generating forces on surface elements. 

The Cauchy stress corresponds to the usual definition of a stress: it gives surface forces 
when contracted on a vector normal to a surface which is prescribed in the physical space. 
It can be shown that the Cauchy stress can be written as a derivative of the energy with 
respect to the deformation gradient tensor: 



nn 1 9W 
T = — 



= V dF 

= F-yl 



(42) 



Unlike the second Piola-Kirchoff stress (i.e. the thermodynamic tension) the Cauchy stress 
does not depend on the choice of a reference state. This equality holds in the limit F_ — > 1: 
in this sense it is an Eulerian derivative of the energy. As we will show below, the second 
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Piola-Kirchoff stress becomes identical to the Cauchy stress in the limit where the current 
and reference configuration are identified. 

In order to relate these two definitions of the stress tensor, we need to be able to switch 
between 7]- and F-derivatives. From the definition of the Green-Saint Venant strain 
tensor, infinitesimal displacements read: 

1 



d 2 2 



(d£ T .£ + £ T .dFj 



(43) 



Since t] is symmetric, in a <i-dimensional problem, it has only d(d + l)/2 independent com- 
ponents, while F_ has d 2 . The rotational degrees of freedom are taken into account by 



considering the antisymmetric infinitesimal displacements: 

1 



: £ 



4ii 



% = - {dE T -E - E T -dg) 



(44) 



which is the generator of infinitesimal rotations, and has d(d—l)/2 independent components. 
From equations and (1I4T). it comes: 



dr) + du = dFj.F 



Using this relation, the differential form of an arbitrary function of F_, A(F_), reads: 



(45) 



dA 



Tr 1 w d ^ T 



Tr (H'( d 2 + d ^)£ _1 

(dA / 



Using the property that the contraction of a symmetric with an antisymmetric tensor van- 
ishes and the fact that r\ is symmetric, we see that the derivative of a function A with respect 
to r] is: 

dA _ 1 
~drj ~ 2 



(46) 



= dg ' dg T ' = 

Incidentally, we also find from the preceding calculation that iff A only depends on t], hence 
does not vary under infinitesimal rotations, it verifies: 



.! dA dA_ T 
'dF dF T '= 



(47) 
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With these formulae in hand, let us come back to the definition (|42|) of the Cauchy stress. 
Taking the limit F_ — > 1 (or t] — > 0), in equation (|4"7)l and looking at equation (|42|) we see that 
T is symmetric. Using equation (j4fij) . it appears that it equals the thermodynamic tension m 
t/iis Hrmt, yet in this limit only. We emphasize here that the symmetry of the stress tensor 
results directly from the invariance of the interaction potential under elementary rotations: 
the energy functional is expected to be a function of {r { } and r\ only: its derivatives with 
respect to the components of u vanish. This is not, however, equivalent to the energy 
functional being invariant under global rotations (the Bravais cell-being a parallelepiped-is 
not). We note, however, that in general dW/dF_ (evaluated off the reference configuration) 
need not be symmetric: this may explain the observation of non-symmetric stresses in 



i on) 

u 



Barron and Klein 



4l| provide the following relation between the Cauchy stress and the 



thermodynamic tension: 

T = —^—F.t.F T . (48) 
= detF = = = v ; 

This relation holds for finite deformations. It can be obtained, following these authors, 

after a Taylor expansion of the energy versus the different strain tensors. We provide here 

a different derivation which relies on an interesting property concerning the transport of 

derivatives. 

As usual, a reference configuration h is given and the system is strained to a current 
configuration h. In order to define the Cauchy stress around the configuration h, we will 
need to consider h as a new reference configuration h = h and also consider a new current 
configuration h . We thus have three different sets of cell coordinates: h,h = h and h . 
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We need to define strain tensors between each pair of these configuration: the deformation 
gradient tensors are denoted F" = hf.h! , Ff = hf.h and F_ = hf.h = h.h . These 
transformations are illustrated schematically in figure Likewise, the Green-Saint Venant 
tensors are denoted: 77", rj and rj. We have the property: Ff = Ff'.F^. For any arbitrary 
function A(F_), and for fixed F_, we can write: 

dA = 




whence, 

dA dA 



).F T . (49) 



dFf' dg 

Using this and equation (J46|) it is then an easy check that: 

™ =F -i™.F-t (50) 
drf = drf = K ' 

This relation simply means that the tensorial derivative with respect to rj transforms as a 
tensor under a change of reference configuration. We just saw in equation (|49jl that this is 
not true of the derivative with respect to F_. 

Going back to the definition of the Cauchy stress (equation (|42j)) and taking the limit 
h — > h = h , we have 

m 1 dW 
T = — 



V OF" 

= F'^F 

and using equation ()50|) in the limit h^h (or F_ — > F), we recover equation 

It is useful to "specialize" the previous expressions for the situations when the energy, or 
any observable can be parameterized as A({rj,- = F_-"tij})- We first have: 

M = (51) 

dg 2 d Lij Uj = 1 ; 

The factor 1/2 results from the fact that each pair is counted twice with the convention of 
implicit summation over repeated indices. Plugging equation (|51|) (and its transpose) into 
equation pUj) we immediately obtain, 

OA 1 , / OA T OA \ T 
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We see from equation ()50|) that this expression consist of a reference-independent formula, 
which is transported backward in the reference configuration as in equation (|5U[). 

Specializing the functional A to A = W, we can now provide an expression for the Cauchy 
stress tensor, using either equation (J48)) and (J52)) or alternatively (|42|) and (J51|) . It reads: 

Z = TT7 ( Sr- z5 + I ti • (53) 



= 4^ V^i~ U 

This is a generalization of the Kirkwood formula to the case of an arbitrary n-body interac- 
tion potential. Note that specialization to systems with pairwise interactions immediately 
yields the standard expression. 

Specializing the function A to A = rfp we can furthermore obtain a useful formula. 
Indeed, in the limit F — > 1, equation (p)2"j) reduces to: 

dr ij 



dr]. 



= 1(8^ + 6^) (54) 



Elastic stiffnesses and constants 

By definition, elastic constants are second derivatives of the energy functional with respect 
to the Green-Saint Venant strain tensor: 

1 d 2 W 



a 



(55) 



V dT] af3 dr] KX 

In this expression, the energy functional W can be any given function of the tensor, rj. For 
example, it can be the energy functional of an atomistic system W(r]) = W({r i:; }, rj) for fixed 
positions of the particles in a reference configuration. Derivatives of such a functional were 
denoted as partial derivatives of U throughout the text. But W could also be defined as the 
energy of an atomistic system following deformation induced changes of a given minimum. 
In this case the partial derivatives of W would correspond to total derivatives for the energy 
functional 1A. The algebra presented here does not depend on these considerations but only 
on the existence of a function W(rj) . 

Since we can commute the order of the derivatives, the elastic constants verify C a p KX = 
Cnxap, and since r\ is symmetric, C a p KX = Cp aKX = C a $ XK . The second order expansion of 
the energy with respect to the components of r] thus reads: 

W - W 1 

— t — = T af3 rj a p + - C aljKX r] aP r] KX + 0{ri a p) (56) 
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where T, 



a3 



■f3a 



(1/V") dW/dF a p\ Q is the stress in the reference configuration. Us- 



ing (|%Tj) and (|56j) . we find for the thermodynamic tension: 

t a /3 = T a j3 + C a/ 3 KX r/ KX + 0(rj 2 x ) 
= T a p + C a/ 3 KX u KX + 0{u 2 KX ) 



(57) 
(58) 



The elastic constant should not be confounded with the quantities which enter the stress- 
strain relations-hence, the wave equation. These are the elastic stiffnesses, c a p KX , and appear 



in the expansion of the energy with respect to tensor u — F_ — l: 



W -W _ o 1 3 

— J-a/3 U a {3 + 77 c af3K X U af3 u k X + V\ U a p) 



(59) 



The relation of the elastic stiffnesses to the elastic constants is obtained by replacing the 
definition (jlljl of r\ in equation (J56j) and expanding in terms of u: 4l| 



c a p KX - C af3KX 2 [2 5 KX T al3 5p x T aK 5 ax 7> K 5p K T ax 5 aK T Px 



(60) 



Note that this expression indicates that elastic stiffnesses do not enjoy the full symmetry of 



4l| 



the elastic constants. Instead, they verify: 

The Cauchy stress can now be obtained from (|57jl by use of equations 



and (HU): |41| 



T a p — T a p + ( c a p KX + T ax 6p K — T a p 5 KX ) Ukx + 0( V 2 kx ) (61) 

Before proceeding to the derivation of an expression for the Born term, let us derive a 
general formula for second derivatives with respect to the components of t]. We first write 
equation (pf6*|) with indices explicitly expressed: 



OA 



dr] 



(r 1 ), 



OA 



+ Or 1 ). 



OA 



y = '«/> dF px v = j xp dF pKA 



The derivative of this expression with respect to T] a a reads: 



d 2 A 



a(£-') 



K,p 



dA 



+ 



» (£-')„ BA 



dr] af3 dF px drj a p dF pK 



+ (r 1 ) 



d 2 A 



*p drj af3 dF px 



we then use the property that in the limit 77 — > 0, 

a (£-'), 



a/3 



F->1 



(r 1 ) 



(62) 



xp dr) aP dF pK 
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by definition, whence we see that 

diF- 1 ) 



after application of equation (j4*fi|. In the limit 77 — »• 0, we thus have for the second derivatives, 
after a further application of equation (jjpj) : 



«9 2 A 



-- 5 



+ - 



&4 



0/3k — 1- 0, 



9 2 A d 2 A 



+ 



+ 



<9 2 A 



+ 



<9 2 A 



(63) 



_9F a pdF KX dF a pdF XK dFg a dF KX dFp a dF XK _ 
We can now provide microscopic expressions for the Born approximation to elastic con- 
stants. We take the energy functional W to be the energy U for some fixed positions of the 
particles in a reference frame. Using equation (J52*j) and (j51j) it is an easy task to write the 
Born term as: 

d 2 u 



/^yBorn 



1 



P X j 

16 V" I ^ ^ K ^ ^ 



d 2 u K d 2 u a x d 2 u 

9rg 0r£ * fci + dri dr!> 13 u 



ij kl 



Q r P Q r X *3 kl 
U 'ij U 'kl 



1 



— j ( T/3 X + S ax Tp K + <5g K T QX + 5p x T a 



(64) 



In order to obtain a similar expression for S iKX , it is convenient to write: 
with 



l-K-X 



d 2 U 



y,«x 



<9rf, 9?7, 



ij "'//ex 



<9 2 W 



77— »0 



17 ' ij 



(65) 



The pair contributions E^. KX are then easily expressed in terms of pair contributions to the 
Hessian: 



1 d 2 U dr 



ki 



2 drZdri dF K 



F->1 



ij w 'kl ~^ K * 

Then, either using equation (JH3|) or applying equation (fH2"J) on /., we find for an 
expression which is analogous to Kirkwood's formula for the stress: 



d 2 W „ K <9 2 W 

^.A _|_ 



ij kl 



ij kl 



(66) 



This general formula permits relating Hj KX to elementary contributions d f Q r of pairs to 
the Hessian. 
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We furthermore notice the similarity between equation (|66|) and equation (|64|). It allows 
writing the following expression for the Born term: 

Ca'pnx = ~~^y ("Sax r ij + S *J>X T ij) ~\ (^ aK ^Px + <W + <W ^X + <W ^) ( 67 ) 

APPENDIX B: MICROSCOPIC EXPRESSION FOR PAIR INTERACTION PO- 
TENTIALS 

We provide here specific expressions for the case when the interaction potential can be 
written as: 

W({e«,-}) = J2 v *i ! ( r «) • 

In this case, the force on individual bonds reads: 

f =- d IlL = J-YlL n .. ( 68 ) 



with, 



r • ■ 

—19 

2tf = 



r 



Introducing the bond tensions and stiffnesses, 



tjj anu Cjj „ , 



the components of the Hessian 

H m 



can be expressed in terms of: 

The elements of the Hessian are then, H = — M for the off-diagonal terms and H 
V ■ , . M for the diagonal terms. 

To obtain an expression for the field S KX , we write: 



E S ^x ( 69 ) 

with (no implicit sum on i and j): 

dr 13 - 1 
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whence, 

"ij'AX = _ ( Ti i Ci 3 ~ n ij n ij U ij ~ 2 ^ ^ ' ft®) 

Inserting this expression in the sum (jHSj) , the second term disappears because the system is 
at mechanical equilibrium. This yields: 

S£* • E ^<j<\r-i,i)»7, »-,»}, ■ (71) 

3 

It is now an easy task to show that the terms in equation (jfi7)) involving contributions 
from the stress tensor disappear, which leads to the expression: 

Cq/3kx = ^ ^ Ci 3 ~ ^ij) r,i 3 U ij U ij U ij ^3 ft"^) 

ij 
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